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To understand the chirality selection in the biological organic system, a simple lattice 
model of chemical reaction with molecular diffusion is proposed and studied by Monte Carlo 
simulations. In addition to a simple stochastic process of conversions between an achiral 
reactant A and chiral products, R and S molecules, a nonlinear autocatalysis is incorporated 
by enhancing the reaction rate from A to R ( or S), when A is surrounded by more than one R 
( or S) molecules. With this nonlinear autocatalysis, the chiral symmetry between R and S is 
shown to be broken in a dense system. In a dilute solution with a chemically inactive solvent, 
molecular diffusion accomplishes chirality selection, provided that the initial concentration 
is higher than the critical value. 
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Organic molecules often have two possible stereo-structures, i.e. a right-handed (R) and 
a mirror-image left-handed (S) form, but those associated with living matter choose only 
one type: only L-amino acids and D-sugards. 1 ' 2 There are many studies on the origin of 
this chirality selection. 1 Various mechanisms proposed to cause asymmetry in the primordial 
molecular environment ( by chance or by external or internal deterministic factors) 2-8 turned 
out to be very minute, and therefore it has to be amplified. 

Frank has shown theoretically that an autocatalytic reaction of a chemical substance with 
an antagonistic process can lead to an amplification of enantiometric excess (ee) and to ho- 
mochirality. 9 Recently, amplification of ee was confirmed in the asymmetric autocatalysis of 
pyrimidyl alkanol, 10-13 and the temporal evolution was explained by the second-order auto- 
catalytic reaction. 12 ' 13 In various other systems such as crystallizations, the chiral symmetry 
breaking is found and discussed extensively. 1 

In our previous paper, 15 we have shown that in addition to the nonlinear autotacalysis a 
recycling process of the reactant introduced by the back reaction accomplishes the complete 
homochirality. There, however, chemical reaction is analyzed macroscopically in terms of aver- 
age concentrations. Thus, a very important factor is neglected, namely the spatial distribution 
of the chemical components; one cannot understand how the homochirality is established over 
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the system. In this paper, we study the chemical reaction of molecules in an extended space 
to understand the proliferation of chirality selection. 
Model and Elementary Processes 

In order to understand the essentials of the chirality selection, we propose here a sim- 
ple model such that the space is restricted to two dimensions and is devided into a lattice. 
Molecules are treated as points moving randomly on the lattice sites. Double occupancy of a 
lattice site is forbidden. There involve four types of molecules in the present minimal model; 
an achiral reactant A, two types of product enantiomers R and S, and a solvant in a diluted 
system. As for the chemical reaction, three typical cases are analyzed; nonautocatalytic, lin- 
early autocatalytic and secondary autocatalytic cases. The back-reaction is always included 
in those three cases. 

The non-autocatalytic chemical reaction proceeds on site and independently of neighbor- 
ing molecules, as illustrated schematically as 



A ^ R 
A ^ S 











R 











The reaction process is essentially stochastic, and we simulate it by the Monte Carlo method. 
In the mean-field approximation where the fluctuation is neglected, the process is described 
by the rate equation in terms of the local concentrations r(i), s(i) and a(i) of R, S, and A 
molecules at a site i as 
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(1) 



with a constant production rate ko and a decay rate A. The rate equation is very powerful in 
theoretical analysis. In this non autocatalytic case, for example, enantiomer concentrations are 
shown to approach to values at a symmetric fixed point = Soo = ^o^oo/A, asymptotically. 
Linearly autocatalytic reaction is described by the reaction scheme 



A + R 
A + S- 
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Since the double occupancy of a lattice site is forbidden in the present lattice model, it is 
natural to assume that these autocatalytic reactions take place when an A molecule is located 
next to one or more R or S molecules, respectively, with a probability k\. Then the additional 
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contribution to the rate equation (1) is described as 

dr(i 
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(2) 



where the summation of sites %\ runs over the 4 nearest neighboring ones to i. 
The nonlinear autocatalysis of the second order is described by the reaction 



A + 2R -> 3R 
A + 2S -> 3S 
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The situation is achieved in our lattice system by assuming that the reaction proceeds when an 
A molecule is surrounded by more than one R or S molecules. By denoting the corresponding 
rate as k 2 , the additional contribution to the rate equation (1) is described as 

dr(i) 

2 



dt 
ds(i) 
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k 2 a(i) ^2 s (*i) s (*2) 



dt 



dr(i) 



dt 



+ 



ds(i) 
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(3) 



where summation over the pairs of sites i\ and %i runs over 6 combinations of the nearest 
neighboring sites to i. 

In addition to the above chemical reaction processes, diffusion of molecules is important if 
the reaction system is diluted in chemically inactive solvent. We assume in the following that 
the diffusion is essentially mediated by the solvent molecules. 
Monte Carlo simulation 

In order to analyze the complex reaction-diffusion system in an spatially extended situa- 
tion, numerical simulation is useful. We adopt a Monte Carlo simulation in two dimensions 
with a following scheme. 

As for the initial condition, A molecules with a concentration cq are distributed randomly 
on a square lattice of a size L x L, and the remaining sites are assumed to be occupied by 
the solvent. Periodic boundary conditions are imposed in the x and y-directions, and the 
length is measured in terms of a lattice constant, hereafter. Then the Monte Carlo simulation 
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starts. One selects randomly a site among L x L square lattice sites. If it is occupied by an A 
molecule, one tries the reaction from A to R with a probability k(A — ► R) and from A to S 
with a probability k(A — » 5). If the A molecule is isolated, k(A — ► R) = k(A — ► 5) = feo- If it 
is surrounded by one R molecules, fc(^4 — ► i?) = fco + fei, and by more than one R molecules, 
k(A — ► R) = ko + k\ + ki- The similar procedure holds when it is surrounded by one or more 
than one S molecules as k(A — > S) = ko + k\ and fc(^4 — > 5) = ko + &i + &2, respectively. If 
the randomly selected site is occupied by an R or S molecule, it is converted to an A molecule 
with a probability A. 

In addition to this reaction algorithm, we can include the diffusion process, if necessary. 
A pair of nearest neighboring sites is chosen randomly, and only if one of them is occupied 
by a solvent molecule, the molecules on the chosen pair sites are exchanged. We assume no 
special bonding among A, R and S molecules, and their diffusion constant is the same, for 
simplicity. Their proximity in space only affects the chemical reaction. Direct exchange among 
A, R and S molecules is excluded so as to realize the diffusionless reaction in the cq = 1 case. 
One Monte Carlo step (MCS) corresponds to L x L trials of chemical reaction and diffusion, 
and thus, on average, each molecule tries reaction and diffusion steps once per each MCS. The 
diffusion constant is D = 1/4 in this time and space unit. 
Results and analysis for the case without diffusion 

For the chirality selection, the macroscopic rate equation shows that the nonlinear auto- 
catalysis plays an essential role. We first perform simulations only with chemical reaction to 
check if our simulation model has the chiral symmetry breaking. The system simulated has a 
size L = 100 with parameters ko = A = 10 -3 , and we change the initial concentration Co and 
the parameters k\ and &2 m order to find out the role of autocatalytic reaction. 

Without diffusion, it is obvious that the result depends strongly on the density Co of the 
reactive molecules. Denser the reactive molecules are, more dominates the nonlinearity. The 
extreme of the dense system is the one without solvent, cq = 1; a whole lattice sites are 
initially occupied by A molecules. At cq = 1 without autocatalysis as k\ = k<i = 0, Monte 
Carlo simulation shows that the numbers of R and S molecules increase synchronously, and 
each saturates about one third of the total lattice sites. The saturation values corresponds 
to the concentrations at the fixed point of the rate equation (1); = Sqo = koCo/(X + 2ko) 
and Oqo = Aco/(A + 2ko), with the initial concentration Co = 1. The asymptotic spatial 
configuration (not shown) is completely irregular, since a molecule on every sites changes its 
state independent of each other. The enantiomatic excess, ee, is defined, as usual, by the 
difference in the concentrations of R and S molecules as 

<P= r — ± - (4) 
r + s 

Without autocatalysis it fluctuates around zero. 

With a linear autocatalysis as k\ = 100/co but = 0, the reaction has produced more R 
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and S molecules( in Fig. 1(b) and (c), respectively), and less A molecule ( in Fig. 1(a)) in a final 
configuration, but the chiral symmetry is not broken; <j> fluctuates around zero, as shown in 
Fig. 2(d). If the spatial dependence is neglected, the corresponding rate equations are written 

as 

dv 

— = (ko + 4&i r)a — Xr 
ds 

— = (k + Akis)a - Xs (5) 

supplimented with the conservation that the sum of the concentrations of A, R and S molecules 
is fixed to its initial value cq; a + r + s = cq. This rate equation has a symmetric fixed point: 
r oo = $oo — (co ~~ a oo)/2 ~ (co/2) — (A/8A4) under the assumption of a strong autocatalysis, 
ko A <C k\. With the present parameters, cq = 1, ko = X = k\/100 = 1CT 3 , the asymptotic 
values = Soo = 0.499 are expected, in fair agreement with the simulation result. Of course, 
in this case, the chiral symmetry should be conserved as the simulation confirms. 

With a second-order autocatalysis with k2 = 100/eo but k\ = 0, the chiral symmetry breaks 
as is shown in Fig.2(a-c). In the simulation shown, there are more R molecules (Fig. 2(b)) than 
S (Fig. 2(c)). By using another sequence of pseudo-random numbers, there occurs equally cases 
that the enantiomer S dominates over R. By neglecting the space dependence the expected 
rate equations are written as 

%- = (k + 6k 2 r 2 )a - Xr 
dt 

ds 

— = (k + 6k 2 s 2 )a - Xs (6) 

with a = cq — r — s. There is a symmetric fixed point U at rjj = su ~ (co/2) — (A/6fc2Co) 
and ajj ~ A/3/C2C0, but it is unstable at a high concentration k 2 c\ S> k^, X. There are 
stable fixed points at Si: (rs 1 ,ss 1 ,as 1 ) ~ (cq — (ko + X)/6k 2 co, k /6k 2 co, X/6k 2 co) and at S2: 
(fS 2 ^ s s 2 i a s 2 ) ~ (ko/Qk 2 co, cq — (ko + X)/6k 2 co, X/6k 2 co). The amplitude of the ee is expected 
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Fig. 2. Configurations of (a) A, (b) R. and (c) S molecules with nonlinear autocatalysis at 10 6 th 
MCS with Co = 1. (d) Evolution of the enantiomeric excess <f>. Upper curve corresponds to the 
case with nonlinear autocatalysis, and the lower curve to that with linear autocatalysis, shown in 
Fig.l. 



to approach to the value given by 

\<f>oo\ ~ 1 - k /3k 2 cl, (7) 

which is close to unity in the present parameter values, = WOko and cq = 1. In the 
simulation, the ee increases gradually as shown by an upper curve in Fig. 2(d), but the final 
asymptotics is not reached yet. 

On looking into the spatial distribution of molecules, one finds that the sites occupied by R 
and S molecules form respective domains. The chirality selection proceeds via the competition 
between two domains. The process is very slow and the configuration in Fig.2(a-c) shows an 
intermediate stage of the S domain shrinking. If both R and S domains extend the whole 
system, the relaxation process becomes even slower. The situation looks quite similar to the 
slow domain dynamics observed in the spinodal decomposition. 

One may notice that the ee in eq.(7) is complete if the nonautocatalytic production is 
absent (ko = 0), in agreement with our previous study. 15 In the present model, we assume 
a finite value of ko, since the creation of an initial chiral molecule R or S from an achiral 
molecule A is necessary. Also without ko, the accidental extinction of the chiral species, R or 
S, is unavoidable and it cannot be recovered only with autocatalytic reactions. 
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(a) 



(b) 



(c) 



Fig. 3. Configurations of (a) A, (b) R, and (c) S molecules with a concentration Co = 0.15 at 3 x 10 4 th 
MCS without diffusion. 



We now consider the reaction and chirality selection in a solution with cq < 1 . The reaction 
system is diluted by adding inactive solvent molecules. If the initial concentration of reactant 
Co is sufficiently high, the nonlinear autocatalysis leads to the chiral symmetry breaking similar 
to the case at cq = 1. On the other hand, at a low concentration far below the percolation 
threshold of the square lattice 16 c p « 0.6, the autocatalysis effect cannot propagate through 
the whole system and fails to break the chiral symmetry. For example, at cq = 0.15, the 
numbers of A, R and S molecules are about the same with each other as shown in Fig. 3(a), 

(b) and (c), respectively. There is no cooperative organization, and the ee fluctuates around 
zero, as shown in Fig. 4(d). 

Results and analysis for the case with diffusion 

Diffusion drastically changes the above situation at low densities. Figures 4(a), (b) and 

(c) , show configurations of A, R and S molecules, respectively, at the concentration cq = 0.15 
with diffusion and nonlinear autocatalysis at the time 3 x 10 4 th MCS. The parameter values 
are D = 1/4, k = A = 10~ 3 , k\ = 0, k 2 = 100fc , and the system size is L 2 = 100 2 . The R 
molecule in Fig. 4(b) increases their number at the cost of A and S molecules in Fig. 4(a) and 
(c). The ee approaches to the saturation value (ftoo ~ 0.73 very quickly. At the fixed point of the 
rate equation (6), the approximate value of the ee, (ft, is given by eq.(7) as 0.871, quite larger 
than the simulation result. However, since k 2 c\ is comparable to ko, the asymptotic form (7) 
is no more valid and the exact value is calculated to be 0.806, closer to the simulation result. 
We have simulated larger system with sizes L 2 = 200 2 and L 2 = 400 2 with k 2 = 100fco, and 
found that the ee becomes non zero as well, though for a large system the initial incubation 
period with vanishing (ft remains very long in some cases. It indicates that there is a certain 
critical size of coherence for the chiral symmetry breaking to take place, and the diffusion 
kinetics controls the propagation of symmetry breaking through the system. 

With a linear autocatalysis and diffusion, our simulation shows that the chiral symmetry 
will not be broken in the diluted system. Therefore, both the nonlinearly autotalytic reaction 
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Fig. 4. Configurations of (a) A, (b) R, and (c) S molecules with a concentration cq = 0.15 at 3 x 10 4 th 
MCS with diffusion, (d) Evolution of the enanthiometric excess (f>. Upper curve corresponds to the 
case with diffusions, and the lower curve to that without diffusion, shown in Fig. 3. 



together with recycling and the diffusion seem to be necessary for the chiral symmetry breaking 
in a dilute solution. 

As we lower the concentration, the value of the ee decreases, and at concentrations lower 
than a critical value about c c ~ 0.12, the system cannot sustain the state with broken chiral 
symmetry. This is also expected from the rate equations (6): it has no symmetry-broken 
solution below the critical concentration 

which is c c = 0.102 for the present choice of the parameters, &2 = 100/co = 100A. 

In fact, this is more clearly understood from the time evolution of the ee, derived from 
the rate equation (6), as 

^ = [ 3 fe 2 ( Co - a) 2 (l - <f) - 2k Q ] (9) 

Here a is the time-dependent concentration of A molecule. The term proportional to &2 i n 
eq. (9) represents that the nonlinear autocatalysis amplifies the chiral symmetry breaking, 
whereas the term proportional to ko suppresses the chiral symmetry breaking by the random 
and independent production of enantiomers, R and S. The state with chiral symmetry looses its 
stability when the coefficient of the linear term in <p in the right-hand side of eq.(9) is positive, 



8/10 



J. Phys. Soc. Jpn. Letter 

and the state with a finite ee, eft, can emerge. Since cq — a represents the total concentration 
of R and S molecules which is close to the initial concentration co, the nonlinear symmetry 
breaking effect becomes weak for a dilute system, and below the critical concentration the 
random creation of racemics dominates. 

As for the critical concentration, there seems to be a discrepancy between the simulation 
result and the rate equation analysis. Since the rate equation corresponds to the mean field 
approximation without fluctuation, the critical concentration in the simulation might turn out 
to be a little larger than the theoretical prediction. Another possibility is the finiteness of the 
diffusion constant. In the rate equation, we assume a homogeneous situation, corresponding 
to an infinitely fast diffusion. With a finite diffusion, the system is influenced by the spatial 
fluctuation. There are also many other possibilities; finite simulation time and size in the Monte 
Carlo simulation, etc. More studies are required on the critical behaviors of this dynamical 
phase transition. 
Summary 

We have proposed a simple lattice model of chemical reaction with molecular diffusion, 
and studied the chirality selection. The nonlinear autocatalysis is shown to be indispensable 
for the selection. In a diluted solution, molecular motion such as diffusion is necessary to 
accomplish the selection. In nature, molecules are in water and the convection should provide 
much more efficient molecular movement. If the initial concentration of the reactant cq is 
too low, below the critical concentration, one can produce only racemic mixture of R and S. 
The critical concentration depends on the ratio of the non-autocatalytic to the nonlinearly 
autocatalytic rate coefficients, k^fki- The asymptotic ee value 4> differs from the complete 
4> = 1 by a factor proportional to k^fki- If the initial production of R or S molecules from 
the reactant A is triggered by minute external effect, fco might be very small and the almost 
complete homochirality be achieved. 

In the whole analysis, the back reaction from the chiral products, R and S, to the reactant 
A is always assumed. Without it, the reaction stops before attaining the full selection, since 
the reactant A is consumed up. The recycling is necessary to develope the selection. But it 
should be smaller than the critical value A c = 2fco (co \/ 6/^2 /ko — 2) in order that the symmetry- 
broken states exist. As A is much too small, the system takes a long relaxation time to reach 
the broken-symmetry state. There seems to be an appropriate range of values of A to achieve 
the chiral symmetry breaking. 
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